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Abstract 

In  this  study,  a  modeling  framework  for  heat  and  mass  transport  is  established  for  a  unit  monolith  type  SOFC,  with  emphasis  on 
quantifying  the  radiation  heat  transfer  effects.  The  Schuster-Schwartzchild  two-flux  approximation  is  used  for  treating  thermal  radiation 
transport  in  the  optically  thin  yttria-stabilized-zirconia  (YSZ)  electrolyte,  and  the  Rosseland  radiative  thermal  conductivity  is  used  to 
account  for  radiation  effects  in  the  optically  thick  Ni-YSZ  and  LSM  electrodes.  The  thermal  radiation  heat  transfer  is  coupled  to  the  overall 
energy  conservation  equations  through  the  divergence  of  the  local  radiative  flux.  Commercially  available  FLUENT™  CFD  software  was 
used  as  a  platform  for  the  global  thermal-fluid  modeling  of  the  SOFC  and  the  radiation  models  were  implemented  through  the  user-defined 
functions.  Results  from  sample  calculations  show  significant  changes  in  the  operating  temperatures  and  parameters  of  the  SOFC  with  the 
inclusion  of  radiation  effects. 
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1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  promising  energy 
source  under  currently  active  development  due  to  their  high 
efficiency,  high  energy  density  and  low  emission  of  pol¬ 
lutants  to  the  environment.  SOFCs  belong  to  the  category 
of  high  temperature  fuel  cells  and  usually  operate  around 
800-1000  °C.  The  high  temperature  operation  allows  for 
the  use  of  natural  gas  as  fuel,  thereby  eliminating  the  need 
for  an  expensive  reformer  system.  Unfortunately,  the  high 
temperature  also  makes  the  fuel  cell  susceptible  to  many 
materials  related  modes  of  thermo-mechanical  failure.  Re¬ 
cently,  mathematical  modeling  has  been  explored  as  a  vital 
tool  in  the  development  and  reliability  prediction  of  fuel 
cells.  An  accurate  modeling  of  the  heat  transfer  mecha¬ 
nisms  in  the  fuel  cell  is  a  prerequisite  to  developing  reliable 
failure  predictions  tools. 

Simplistic  thermal  models  of  fuel  cells,  using  heat  trans¬ 
fer  correlations  for  assumed  fluid  flow  configurations,  have 
been  carried  out  in  the  past  [1].  More  recently,  temperature 
distributions  in  three-dimensional  models  of  planar  SOFC 
have  been  achieved  by  coupling  the  electrochemical  effects 
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with  the  solution  of  the  Navier-Stokes  equations  of  motion 
and  energy/mass  conservation  [2]. 

Most  existing  models  neglect  the  effect  of  radiation  heat 
transfer.  Yakabe  et  al.  [3]  found  that  inclusion  of  radia¬ 
tion  heat  exchange  within  the  flow  channels  resulted  in 
a  flatter  distribution  of  the  temperature  profile  along  the 
fuel  cell.  However,  the  above  model  neglects  absorption, 
emission  or  scattering  in  the  media  and  only  accounts 
for  the  surface-to-surface  radiation  effects  in  the  analysis. 
Though,  Hartvigsen  et  al.  [4]  were  able  to  demonstrate 
significant  changes  in  stack  temperature  with  the  inclusion 
of  radiation  effects,  not  much  detail  was  provided  in  terms 
of  the  properties  used  and  the  specific  radiation  model 
employed. 

In  the  present  study,  the  effect  of  radiative  heat  transfer  on 
the  fuel  cell  performance  is  reported.  A  single  unit,  mono¬ 
lith  type,  co-flow  fuel  cell  (Fig.  1)  model  was  studied  by  an¬ 
alyzing  the  fluid  flow,  species  transport  and  conjugate  heat 
transfer  within  the  cell.  The  radiative  heat  transfer  (RTE) 
equation  was  used  to  model  radiation  heat  transfer  and  was 
solved  using  the  discrete  ordinate  (DO)  method.  Alternately, 
a  simplified  radiation  model  was  also  developed  based  on 
the  Schuster-Schwartzchild  two-flux  approximation  and  its 
accuracy  was  compared  to  that  of  the  computationally  ex¬ 
pensive  DO  scheme. 
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Nomenclature 

c o  speed  of  light  in  vacuum  (m/s) 

G  Gibb’s  free  energy  (J) 

h  Planck’s  constant  (J  s) 

H  reaction  enthalpy  (J) 

k  Boltzmann  constant  (J/K) 

k r  Rosseland  thermal  conductivity  (W/(mK)) 

n  refractive  index  of  the  medium 

<7r  radiative  heat  flux  (W/m2) 

T  absolute  temperature  (K) 

Tr  transmittance 

Greek  letters 

/Jr  Rosseland-mean  extinction  coefficient  (m-1) 
k  absorption  coefficient  (m-1) 

X  wavelength  (|xm) 

a  Stefan-Boltzmann  constant  (J/(K4  m2  s)) 

r  optical  thickness 


Cathode 


2.  Modeling  methodology 

2.1.  Radiative  properties 

The  prediction  of  radiation  heat  transfer  in  semi-trans¬ 
parent  materials  requires  the  knowledge  of  radiative  prop¬ 
erties,  namely,  the  absorption  coefficient  and  the  scattering 
coefficient  of  the  material.  The  state-of-the-art  SOFCs 
uses  yttria-stabilized-zirconia  (YSZ)  as  the  electrolyte, 
strontium-doped  LaMnC>3  (LSM)  as  the  cathode  and  nickel 
doped  YSZ  as  the  anode.  The  radiative  properties  of  a  ma¬ 
terial  are  usually  a  function  of  wavelength,  and  the  spectral 
region  of  interest  for  thermal  radiation  can  be  obtained  from 
the  Planck’s  law  for  blackbody  emissive  power, 

2nhcQ 

Ebk  =  n2X5[ehc°/nkkr  -  1]  (1) 

where  is  the  spectral  emissive  power  of  a  black  body,  h 
the  Planck  constant,  co  the  speed  of  light  in  vacuum,  k  the 
Boltzmann’s  constant,  T  the  temperature  and  n  the  refrac¬ 
tive  index  of  the  medium.  The  variation  of  the  black  body 


Fig.  2.  Spectral  black  body  emissive  power  for  n  =  1.6,  T  =  973  K,  80% 
of  radiative  energy  contained  within  1.4  p,m  <  /.  <6.1  p.m. 


emissive  power  with  wavelength  for  n  =  1.8  [5]  and  T  = 
973  K  is  shown  in  Fig.  2  and  corresponds  to  SOFC  operating 
conditions.  It  can  be  shown  that  80%  of  the  radiative  emis¬ 
sive  power  is  contained  within  the  infrared  spectral  region 
1.4  |xm  <  X  <  6.1  p,m  with  the  maximum  at  X  ^  1.86  p,m. 

Fig.  3  shows  the  radiative  properties  of  YSZ  obtained  with 
a  500  p,m  thick  sample.  Clearly,  a  gray  body  approximation 
can  be  used  in  the  radiation  modeling  of  YSZ,  as  the  trans¬ 
mittance  is  found  to  be  independent  of  wavelength  in  the 
region  of  interest  (1.4  p,m  <  X  <6.1  p,m).  The  absorption 
coefficient  of  YSZ  electrolyte  is  obtained  from 

Tr  =  q~kL  (2) 

where  Tr  is  the  transmitivity,  k  the  absorption  coefficient  and 
L  the  thickness  of  the  medium.  The  optical  properties  of  the 
cathode  are  obtained  from  the  imaginary  part  of  dielectric 
constant  spectra  [6]  and  the  optical  conductivity  spectra  [7] 
of  Lai_  vSrAMn03. 

The  optical  thickness  (r  =  kL)  of  the  fuel  cell  electrolyte 
and  cathode,  based  on  the  dimensions  in  Table  1,  were  found 


Fig.  3.  Optical  properties  of  a  500  |xm  thick  YSZ  sample. 


S.  Murthy,  A.G.  Fedorov  /  Journal  of  Power  Sources  124  (2003)  453—458 


455 


Table  1 

Monolith  fuel  cell  model  description 


Region 

Type 

Dimension  (m) 

Cathode 

Porous  (0 

=  0.3) 

0.0025 

Anode 

Porous  (0 

=  0.38) 

0.0025 

Electrolyte 

Solid 

0.0005 

Current  collector 

Solid 

0.0005 

to  be  0.25  and  10  ,  respectively.  The  optical  properties  of  the 
anode  material  are  complicated  to  calculate  as  they  depend 
on  the  distribution  of  the  nickel  cermets  in  YSZ.  The  pres¬ 
ence  of  Nickel  doping  makes  the  anode  highly  absorbing  to 
thermal  radiation,  thus  resulting  in  optically  thick  (r  »  1) 
behavior.  Therefore,  for  lack  of  experimental  data,  the  ab¬ 
sorption  coefficient  of  the  anode  is  assumed  to  be  the  same 
as  that  of  the  cathode.  Also,  the  radiation  scattering  in  elec¬ 
trodes  is  neglected,  although  this  assumption  is  questionable 
considering  the  porous  nature  of  electrodes. 

2.2.  Radiative  model  approximation 


2.2.1.  Rosseland  thermal  conductivity 

For  optically  thick  slabs  (r  1),  it  has  been  shown  [8] 
that  the  radiative  heat  flux  q r  can  be  estimated  with  sufficient 
accuracy  through  Rosseland  approximation  as, 

qR  =  -fRvr  (3) 


where  kR  is  the  radiative  conductivity  and  defined  as, 


kR  = 


16«2aT3 

3  At 


(4) 


and  lower  hemisphere.  For  a  gray,  non-scattering  medium 
confined  between  two  isothermal,  parallel  black  plates  at 
temperatures  7’top  and  '/’bottom  anc*  separated  by  a  distance 
L,  the  two-flux  model  gives  the  radiative  heat  flux  as  [8] 

qR(z)  =  C  i  elKZ  +  C2e~2KZ  (5) 

where 

Cl  =  -a(T4p  -  r4)e-2^ ,  C2  =  a(rb4ottom  -  T4)  (6) 

The  radiation  transport  thus  calculated  can  be  coupled 
with  the  overall  energy  conservation  by  introducing  the  di¬ 
vergence  of  the  radiative  heat  flux  as  a  (negative)  source 
term  in  the  energy  equation. 

Fig.  4  shows  a  schematic  of  a  two-dimensional  enclosure 
used  for  validating  the  Schuster-Schwartzchild  radiation 
model.  The  bottom  and  top  walls  of  the  enclosure  were 
maintained  at  constant  temperatures  of  1000  and  2000  K, 
respectively,  with  the  sidewalls  being  adiabatic.  Studies 
were  carried  out  for  different  optical  thickness  by  varying 
the  absorption  coefficient  of  the  enclosed  medium. 

Fig.  5  shows  the  temperature  variation  across  the  thick¬ 
ness  of  the  liquid  medium  as  predicted  by  the  two-flux 
approximation,  and  compares  their  accuracy  with  those 
obtained  using  the  more  accurate  but  computationally  in¬ 
tensive  DO  scheme.  Excellent  agreement  in  temperature 
predictions  is  seen  for  extremely  optically  thin  materi¬ 
als.  Large  deviations  are  observed  with  increase  in  optical 
thickness  with  the  Schuster-Schwartzchild  approximation 
failing  for  optical  thickness  greater  than  one. 

2.3.  Fuel  cell  model 


In  Eq.  (4),  /3R  refers  to  the  Rosseland-mean  extinction  coef¬ 
ficient,  which  is  equal  to  the  constant  absorption  coefficient 
for  the  gray  and  non-scattering  medium.  Unfortunately,  the 
above  diffusion  approximation  is  not  valid  near  a  bound¬ 
ary,  where  it  often  fails  quiet  miserably  [8],  This  limitation 
can  be  overcome  by  coupling  the  diffusion  model  with  the 
Schuster-Schwartzchild  two-flux  approximation  for  locally 
optically  thin  regions  near  the  boundaries  of  the  overall  op¬ 
tically  thick  medium. 

2.2.2.  Schuster-Schwartzchild  two-flux  approximation 
Schuster-Schwartzchild  or  the  two-flux  approximation 
provides  a  simple  solution  method  for  one-dimensional 
radiation  in  a  thin,  plane-parallel  slab.  In  this  method,  the 
radiative  intensity  is  assumed  to  be  an  isotropic  function 
of  the  propagation  direction,  but  different,  over  the  upper 


A  general-purpose  computational  fluid  dynamics  soft¬ 
ware  (FLUENT™)  was  used  for  modeling  the  transport 
and  heat  transfer  phenomena  occurring  within  the  fuel  cell 
[9],  The  Navier-Stokes  and  transport  equations  are  solved 
to  obtain  the  flow,  species  concentrations  and  temperature 
at  each  location  in  the  cell.  A  coupled  electrochemical 
model  calculates  the  current  density,  cell  voltage,  and  heat 
flux  at  the  electrodes  based  on  the  species  and  temperature 
distribution  within  the  cell.  This  information  is  then  passed 
back  to  CFD  model,  where  it  is  used  to  update  the  species 
and  temperature  fields.  The  process  continues  in  an  iterative 
manner  until  convergence  is  reached. 

The  single-cell  model  (Fig.  1)  represents  a  repeating  unit 
cell  in  the  center  of  a  large  stack.  Hence,  adiabatic  bound¬ 
ary  conditions  were  assumed  on  all  cell  walls.  Four  chem¬ 
ical  species,  namely  H2,  02,  H20  and  N2,  were  tracked  in 


r  Tjjp  =  2000  K 


^Centerline  -  temperature  comparison 


Non-scattering  medium 

A - - — - ► 

1  cm 


Tbcttora  —  1000  K " 


20  cm 


Fig.  4.  Schematic  of  two-dimensional  model  in  FLUENT™  used  for  validation  of  Schuster-Schwartzchild  two-flux  radiation  model. 
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Fig.  5.  Comparison  of  two-flux  method  and  more  accurate  discrete  ordinate  (DO)  scheme  for  different  optical  thickness.  The  two-flux  method  fails  for 
optical  thickness  greater  than  unity. 


the  simulation.  The  inlet  flows  of  air  and  fuel  were  speci¬ 
fied  using  constant  mass  flow  boundary  conditions.  Air  was 
delivered  to  the  cathode  at  0.32  mg/s  and  750  K.  The  fuel 
mass  flow  rate  was  0.0045  mg/s  at  750  K  with  the  composi¬ 
tion  of  fuel  by  mass  fraction  being  78%  kb  and  28%  H2O. 
A  fixed  pressure  condition  was  used  for  the  fuel  and  air  out¬ 
lets.  The  pressure  drops  in  the  porous  cathode  and  anode 
regions  were  modeled  using  Darcy’s  law,  and  the  transport 
of  the  reactant  species  to  the  electrolyte/anode  interface  was 
modeled  through  diffusion  of  a  multi-component  mixture 
through  porous  electrodes. 

The  air  and  fuel  streams  were  modeled  as  ideal  gas  mix¬ 
tures  with  the  gas  viscosities  and  thermal  conductivities  cal¬ 
culated  based  on  kinetic  theory.  The  material  properties  of 
the  solids  used  in  the  calculations  are  listed  in  Table  2.  The 
anode  and  cathode  regions  were  modeled  as  porous  media 
with  porosities  of  0.38  and  0.30,  respectively.  The  material 
properties  of  the  porous  medium  were  computed  as  the  vol- 


Table  2 

Solid  material  properties 


Material 

Density  (kg/m3) 

Thermal  conductivity 
(W/(m  K)) 

Anode-solid 

3030 

5.84 

Cathode-solid 

3310 

1.86 

Electrolyte 

5160 

2.16 

Current  collector 

8030 

0.1 

ume  average  of  the  fluid  and  solid  properties.  The  three  di¬ 
mensional  simulations  of  the  unit  cell,  including  the  inter¬ 
nal  gas  flow  channels,  were  performed  using  a  fine  mesh  of 
37,440  computational  elements. 


3.  Results 

3.1.  Radiation  effect 

Fig.  6  shows  the  temperature  variation  along  the 
cathode-electrolyte  and  anode-electrolyte  walls  of  the  fuel 
cell  with  and  without  inclusion  of  radiation  effects.  Clearly, 
the  heat  transfer  enhancement  due  to  radiation  results  in 
overall  decrease  of  the  temperature  level  within  the  cell,  as 
well  as  the  less  steep  temperature  gradients  along  the  cell. 

The  maximum  electrical  work  obtainable  in  a  fuel  cell, 
operating  at  constant  temperature  and  pressure,  is  given  by 
the  change  in  Gibb’s  free  energy  (AG)  of  the  electrochem¬ 
ical  reaction.  The  relationship  between  the  Gibb’s  free  en¬ 
ergy  (G)  and  the  reaction  enthalpy  (H)  is  known  to  be  [10], 

AG  =  AH  —  T  AS  (7) 

where  A H  corresponds  to  the  total  thermal  energy  available 
in  the  system,  AS  denotes  the  change  in  entropy  and  T  AS 
represents  the  amount  of  heat  generated  by  a  fuel  cell  oper¬ 
ating  reversibly.  The  change  in  entropy  (AS)  is  a  non-linear 
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Fig.  6.  Temperature  along  a  centerline  passing  through  the  length  of  the 
fuel  cell  (^-direction)  on  (a)  cathode-electrolyte  and  (b)  anode-electrolyte 
walls  with  and  without  accounting  for  radiation  heat  transfer. 

function  of  temperature  and  is  of  the  form, 

A  S=  AS°  +  aln^^J  +  b(T  -  298)  +  ^c(T  -  298)2 

(8) 

where  a,  b ,  and  c  are  empirical  constants  [10]  and  A.S’°  is 
the  standard  entropy  of  formation  at  298  K. 

The  additional  heat  transfer  path  provided  by  radiation 
results  in  a  fairly  uniform  temperature  field  within  the  fuel 
cell  system  with  reduction  in  the  maximum  temperatures. 
Inclusion  of  radiation  effects  resulted  in  a  20%  drop  in  the 
T  AS  value,  thereby,  accounting  for  the  1 80  K  decrease  in  the 
fuel  cell  temperature.  Also,  the  maximum  available  electrical 
work,  as  given  in  Eq.  (7),  increases  with  decrease  in  T  AS 
and  results  in  corresponding  higher  fuel  cell  voltage.  Indeed, 
the  observed  180  K  drop  in  the  cell  temperature  (Fig.  6) 
causes  a  corresponding  increase  in  the  Gibbs  free  energy, 
and  results  in  a  significant  increase  in  the  cell  voltage  from 
0.65  to  0.74  V. 

3.2.  Rosseland/two-flux  model 

Although  the  DO  method  is  very  accurate  in  accounting 
for  absorption  and  emission  of  thermal  radiation  through  the 


solution  of  the  RTE,  the  computational  costs  involved  along 
with  the  memory  requirements  make  the  implementation 
prohibitive  for  extensive  design  calculations  and  optimiza¬ 
tion.  To  overcome  the  above  problem,  a  simplified  radiation 
model  has  been  developed  that  takes  advantage  of  the  dif¬ 
ferences  in  optical  properties  of  the  fuel  cell  materials  to  re¬ 
duce  the  computational  expense.  In  this  approach,  radiative 
heat  transfer  in  the  optically  thin  electrolyte  was  modeled 
using  the  Schuster-Schwartzchild  model,  whereas  radiation 
in  the  optically  thick  electrodes  was  accounted  through  the 
Rosseland  diffusion  approximation.  Failure  of  the  diffusion 
model  near  the  boundaries  was  overcome  by  employing  a 
very  fine  grid  near  the  boundaries.  Subsequently,  the  first  few 
cells  along  the  boundary  walls  of  the  electrode  were  treated 
as  optically  thin  medium  and  solved  using  the  two-flux  ap¬ 
proximation  method  as  well. 

Fig.  7  shows  the  variation  of  temperature,  along  the  length 
of  the  anode-electrolyte  interface  (x-direction)  and  across  the 
thickness  of  the  fuel  cell  (at  the  location  x  —  0.025  m),  for 
the  different  radiation  models.  The  maximum  difference  in 
the  temperature  predictions  obtained  using  a  highly  accurate 
discrete  ordinate  (DO)  method  was  found  to  be  <10K  as 
compared  to  the  simplified  Rosseland/two-flux  approximate 
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Fig.  7.  Validation  of  Rosseland/two-flux  radiation  sub-model  through 
temperature  comparison  with  DO  predictions  along  (a)  centerline  passing 
along  the  anode-electrolyte  interface  and  (b)  thickness  of  the  fuel  cell  (at 
axial  position  x  =  0.025  m). 
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scheme.  At  the  same  time,  the  latter  radiation  sub-model 
resulted  in  almost  a  ten-fold  reduction  of  the  CPU  time 
from  658  min  for  the  DO  method  to  only  76  min  with  the 
simplified  Rosseland/two-flux  approximation  scheme. 

4.  Conclusions 

Based  on  the  rigorous  analysis  of  radiation  heat  transfer  in 
a  single  unit  SOFC,  the  following  conclusions  can  be  drawn: 

(i)  Though  commonly  neglected,  radiation  heat  transfer 
effects  are  significant  and  need  to  be  accounted  for 
accurate  prediction  of  the  temperature  field  and  the  fuel 
cell  output  voltage. 

(ii)  The  computationally  expensive  discrete  ordinate 
method  can  be  replaced  by  a  simplified  radiation 
sub-model,  which  is  based  on  Rosseland  diffusion 
approximation  for  the  optically  thick  electrodes  and 
the  Schuster-Schwartzchild  two-flux  method  for  the 
optically  thin  electrolyte. 

(iii)  Significant  saving  in  computational  time  is  achieved 
with  the  Rosseland/two-flux  radiation  model  without 
much  loss  in  accuracy  in  the  temperature/voltage  pre¬ 
dictions  when  compared  to  the  discrete  ordinate  (DO) 
method. 
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